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The non-trivialness of a topological insulator (TI) is characterized either by a bulk topological 
invariant or by the existence of a protected metallic surface state. Yet, in realistic samples of 
finite size this non-trivialness does not necessarily guarantee the gaplessness of the surface state. 
Depending on the geometry and on the topological indices, a finite-size energy gap of different nature 
can appear, and correspondingly, exhibits various scaling behaviors of the gap. The spin-to-surface 
locking provides one of such gap-opening mechanisms, resulting in a power-law scaling of the energy 
gap. Weak and strong TI's show different degrees of sensitivity to the geometry of the sample. As a 
noteworthy example, a strong TI nanowire of a rectangular prism shape is shown to be more gapped 
than that of a weak TI of precisely the same geometry. 



I. INTRODUCTION 

The non-trivialness of a TI is often characterized by the 
presence of a gapless surface state. 1 ' 2 A one-to-one corre- 
spondence can be established between the (non-) trivial- 
ness of a bulk topological invariant and the presence vs. 
absence of the gapless surface state (bulk-surface corre- 
spondence) . However, precisely speaking, for such a gap- 
less state to be existent, both the trivial and non-trivial 
sides are semi-infinite, separated by an infinitely large 
interface. The above distinction can be made, therefore, 
only in such an idealized situation. TI samples, in real- 
ity, occupy only a finite domain of the space, and have 
also a variety of shapes surrounded generally by a curved 
or folded surface(s). In experiments it is also the case 
that some TI samples of nanometer scale size exhibit a 
clear gapless surface state, while other samples of the 
same chemical composition but of a different geometrical 
shape do not necessarily exhibit a clear signature of topo- 
logical non-triviality. Such an issue will be addressed in 
this paper. 

The main scope of the paper is concomitant with the 
observation that there are three different gap-opening 
mechanisms effective in the samples of finite size. The 
most primitive among them is the one due to mixing 
of the surface electronic wave functions on the opposing 
sides, e.g., of an infinitely large slab-shaped sample. Such 
an energy gap associated with the finite thickness of the 
gapped bulk, decays exponentially as a function of the 
thickness of the slab, and is in practice almost irrelevant 
except in extremely thin film samples. 3 The low-energy 
(surface) electronic spectrum in the slab geometry suffers 
indeed only from this type of exponentially small finite- 
size energy gap. 4-7 The second mechanism to open a gap 
in the surface electronic spectrum, which is also more 
relevant in magnitude, is the so-called "spin-to-surface 
locking" . 8-11 The electronic spin in the a priori gapless 
surface state on a curved surface of TI has a tendency 
to be locked in-plane to the local tangent of the surface. 
In the cylindrically symmetric case, the spin-to-surfacc 
locking results in the half-integral quantization of the or- 



bital angular momentum along the axis of the cylinder. 
The half-odd integer quantization gaps out the spectrum, 
and this gap decays only algebraically; qualitatively more 
relevant than the gap of the previous type. The spin-to- 
surfacc locking leads, indeed, irrespective of the presence 
of cylindrical symmetry, e.g., in a prism-shaped sample, 
to opening of the gap. 

Another aspect of the topological insulator which we 
aim at exploring in this paper is the role of anisotropy, 
especially in the weak topological insulator (WTI) phase. 
This is much related to the third mechanism of gap- 
opening, which occurs due to the interplay of the 
anisotropy of WTI and the specific geometry we will fo- 
cus on (case of the prism-shaped geometry). In three spa- 
tial dimensions (3D), Z 2 topological insulator is known 
to be characterized by four Z2 indices, 12-14 the principal 
(strong) index v and other "weak" indices ^1,^2,^3, in- 
stead of a single Z 2 index in the case of 2D. The principal 
index vq is used to distinguish a STI (y§ = 1) from trivial 
and weak topological insulators (v = 0). In a WTI at 
least one of its weak indices exhibits a nontrivial value 
(=1). A WTI shows generally an even number of helical 
Dirac cones on its surfaces, but on the surface normal 
to its "weak vector" V = (z^i, ^2, ^3) it shows no Dirac 
cone. The WTI can be viewed as stacked layers of 2D Z2 
topological insulators. In this regard the set of weak in- 
dices (y\Viv-i) can be regarded as the Miller index of such 
stacked layers. Since gapless surface states are expected 
to form only at the edge of the stacked layers, one can 
naturally understand that no Dirac cone is formed on a 
surface normal to v in this picture. To summarize, the 
WTI bears two Dirac cones on surfaces parallel to V and 
no Dirac cone on surfaces normal to V. When this char- 
acteristic feature is combined with the specific (rectangu- 
lar) prism geometry, the anisotropy of a WTI manifests 
as an alternating size-dependence of the energy gap; the 
magnitude of the gap is qualitatively different whether 
the number of "stacked layers" is even or odd. It will be 
demonstrated that weak and strong TI's show different 
degrees of sensitivity to the geometry of the sample. 

The periodic table of topological insulators and super- 
conductors classifies them by the nature of strong in- 
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dices characterizing the system. The weak indices are 
not shown, at least explicitly on the table. 15,16 Showing 
an even number of Dirac cones on its surfaces, WTI is a 
priori considered not to be robust. But recently, a few 
counter examples to this common belief have been pro- 
posed. One is the existence of protected gapless helical 
modes along a dislocation line in the WTI. 17 ' 18 More re- 
cently, a couple of papers have appeared, demonstrating 
that an even number of Dirac cones on the surface of WTI 
are actually not that fragile against disorder. 19,20 Here, 
we point out that in a specific situation in the prism- 
shaped geometry the surface state of a WTI is in a sense 
"more strongly protected" from a finite-size energy gap 
than that of a STL 

The paper is organized as follows. In Sec. II we intro- 
duce our effective model Hamiltonian for 3D anisotropic 
topological insulators. The phase diagram of the model 
is determined by the calculation of topological numbers 
in the bulk. In Sec. Ill we discuss different origins of 
the finite-size energy gap, highlighting the role of spin- 
to-surface locking in the cylindrical geometry. In Sec. VI 
we demonstrate that in the more realistic rectangular- 
prism geometry, three types of gap-opening appear and 
disappear by a small change of model parameters, lead- 
ing to an intricate size dependence of the gap. Sec. VI 
is devoted to Conclusions. 



II. MODEL AND ITS PHASE 
DIAGRAM— ENGINEERING THE WEAK 
INDICES 

As a concrete realization of strong and weak topo- 
logical insulators with specific strong and weak indices, 
Vq and v = (1/1,1/2,^3), we consider as given in Eq. 
(1), a Wilson-Dirac type effective Hamiltonian for a 3D 
topological insulator implemented on a cubic lattice. 21,22 
Since we will be interested in the analysis of WTI phases 
with anisotropic weak indices, we choose the mass pa- 
rameters 77122;, rn,2y, rri2z appearing in the Wilson term 
[see Eq. (2)] to be anisotropic. 

A. The Wilson-Dirac type effective Hamiltonian 

Let us consider the following Wilson-Dirac type effec- 
tive Hamiltonian for a 3D topological insulator imple- 
mented on a cubic lattice: 

-Hbuik = e(fc)l + T x m(k) + TyCr^A^ sinfc M , (1) 

where e(fc) is an even function of fc, and 

m(k) = too + 2m 2A1 (l - cos fc M ). (2) 

In Eqs. (1) and (2) a summation over the repeated index 
/i (= x, y, z) is not shown explicitly. The model specified 
by this couple of equations can be regarded as a tight- 
binding model with only the nearest neighbor hopping, 
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FIG. 1. The phase diagram of the Wilson-Dirac type effective 
tight-binding Hamiltonian given in Eqs. (1), (2). Notice the 
anisotropy of our hopping parameters [see Eqs. (3)]. In each 
of the strong (STI) and weak (WTI) topological insulator 
phases, together with the nature of the specific phase, the 
four Z2 indices Vj (j = 0,1,2,3) and the winding number 
N3 are shown, as Nz{vo,v\V2Vd,). The solid lines representing 
the phase boundaries correspond to closing of the bulk energy 
gap. 

determining the structure of the energy bands over the 
entire 3D Brillouin zone (BZ). Eq. (1) can be regarded 
as a 4 x 4 matrix, spanned by two types of Pauli ma- 
trices a and r each representing physically real and or- 
bital spins, respectively. Compared with a more generic 
representation of the Dirac Hamiltonian in terms of the 
"7-matrices" , we have chosen in Eq. (1) "70" coupled to 
the mass term m{k) associated only with an orbital spin 
r x . 

The mass term (2) represents (a half of) the band gap 
at time-reversal invariant momenta (TRIM), k = fco, sat- 
isfying — ko = fco + G, with G being a reciprocal lattice 
vector, corresponding either to a normal or an inverted 
gap, depending on the relative sign of Too and the coef- 
ficient of the quadratic (Wilson) term at a given TRIM. 
By investigating this feature of band inversion at the 
eight TRIM as varying the mass parameters, one can 
identify 23 various weak and strong TI phases character- 
ized by strong and weak indices, vq and v = (1/1, V2, 1/3). 
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Phase boundaries between such topologically distinguish- 
able insulating phases correspond necessarily to closing 
of the bulk energy gap. 

Known examples of 3D topological insulators are lay- 
ered materials, exhibiting, in the leading order approxi- 
mation, uniaxial anisotropy in the crystal c-axis. 24-28 To 
reflect this feature in the effective tight-binding model, 
i.e., in Eqs. (f), (2) we assume that our model parame- 
ters have the same uniaxial anisotropy. 21 Especially, the 
three mass parameters ni2x, m2y, ~m 2z are classified to 
two types: m 2 y and m 2 j_, depending on whether the cor- 
responding hopping direction is, either parallel or per- 
pendicular to the stacked layers of the crystal. Clearly, 
the correspondence depends on the relative orientation 
of the crystal growth axis and our cartesian coordinates; 
e.g., when the crystal c-axis is oriented to the direction 
of z-axis, 

m 2J _ = m 2z , m 2 \\ = m 2x = m 2y , 
A ± = A Z , A,, = A x = A y . (3) 

Independently of this choice of the relative orienta- 
tion, our control parameters for specifying topologically 
different phases are relative magnitudes of mo, m 2 ± 
and m 2 i|. Then, by studying the feature of band in- 
version at eight TRIM as a function of these control 
parameters, 23 one can deduce the phase diagram of the 
model. FIG. 1 shows such a phase diagram depicted in 
the (m /m 2 ||,m 2 _L/m 2 ||)-plane. 



B. Phase diagram 



Stopping at mo/ m 2|| — — 1, let us now vary m 2 j_/m 2 \\ , 
i.e., introduce anisotropy in the mass parameters. On 
the line m^/m 2 \\ = —1 (a thick red line in FIG. 1), the 
system is in a STI phase with i/q = 1 and V = (0,0,0), 
when TO 2 j_/m 2 || > 1/4. The anisotropy appears in the 
weak indices below this critical value, m 2 ± = — mo/4, 
corresponding to band crossing occurs at the Z-point, 
and the system enters a WTI-A phase with v = and 
v = (0,0,1) when m 2 j_/m, 2 || < 1/4. In later sections 
we will quantify various manifestations of this quantum 
phase transition in the finite size effects. The situation 
is similar on the line mo/m 2 j| = —5 (a thick blue line in 
FIG. 1), above and below the critical point m 2 j_/m, 2 || = 
1/4, although in this second example the transition oc- 
curs from an isotropic to an anisotropic WTI phase, each 
named, respectively, WTI-B and WTI-C phases. 
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TABLE I. Definition of the surfaceless, slab, (rectangular) 
prism and cubic geometries. Here, to avoid confusion in the 
terminology, we define these different types of geometries in 
terms of the switching on and off of the periodic boundary 
conditions (PBC) in the x-, y- and z-directions. In the Ta- 
ble, "1" and "0" signifies that the PBC in the corresponding 
direction is, respectively, on and off. In the latter case, PBC 
is replaced by the fixed boundary condition (FBC). 



FIG. 1 shows the phase diagram of the Wilson-Dirac 
type effective tight-binding Hamiltonian given in Eqs. 
(1), (2). The uniaxial anisotropy of the hopping param- 
eters, as given by Eqs. (3), is taken into account. Each 
of the STI and WTI phases are characterized by four 
Z 2 indices. The calculated winding number N 3 (see Ap- 
pendix A) is also shown. Solid lines, separating neigh- 
boring topologically distinct phases, indicate closing of 
the bulk energy gap. Duplicate lines appearing at the 
phase boundary correspond to simultaneous formation of 
two bulk 3D Dirac cones. The duplication is due to the 
uniaxial choice of the hopping parameters. To see such 
specific features, let us focus below on a few particular 
examples of the STI and WTI phases. 

Let us first concentrate on the isotropic line 
m 2J_/ m 2\\ = 1 in the phase diagram (indicated as a thick 
green line in FIG. 1). The change of the winding number 
A*3 on this line is shown in the first panel of FIG. 11. No- 
tice that on this line different STI and WTI phases show 
only symmetric weak indices. At the phase boundaries 
between STI and WTI phases, a double and single solid 
lines cross, indicating simultaneous closing of three Dirac 
cones in the bulk. This occurs at A, Y, Z: kx = (it, 0, 0), 
k Y = (0,7r,0), k z = (0,0, 7r), three symmetric points 
(TRIM) in the 3D BZ. 



III. DIFFERENT ORIGINS OF THE 
FINITE-SIZE ENERGY GAP 

A single Dirac cone on the surface of a STI is topolog- 
ically protected, 2 and also robust against disorder. 29 ' 30 
In reality, TI samples always have a finite thickness be- 
tween the two surfaces of opposing sides. Imagine a slab- 
shaped sample (c./. Table I), which we assume infinitely 
large, neglecting the existence of side surfaces. In such a 
slab geometry, STI bears a pair of surface Dirac cones, 
each localized in the vicinity of the two opposing surfaces. 
These two "Dirac cones" do not communicate, and con- 
sequently remain gapless, as far as the thickness of the 
slab is much larger than the penetration of the surface 
state into the bulk [see Appendix B for an extensive dis- 
cussion on the penetration of the surface wave function 
in the slab geometry; see also Refs. 4 ~ 7 ]. 

In a sense this gaplessness is also protected by the very 
slab geometry. In the case of a sample of more realistic 
shape with typically side surfaces (c/., cases of a prism 
and a cube; see Table I), the same protection is no longer 
valid. The side surfaces open a priori gapless channels 
allowing for communication between the two initial Dirac 



4 



cones on two surfaces of the slab. Since this communica- 
tion through gapless side surfaces is much stronger than 
the one through the gapped bulk (cf. case of the slab 
geometry) , it leads to opening of a size gap qualitatively 
more relevant than the latter case. 

Of course, the effects of such side surfaces appear in 
the transport characteritics only when an electron can 
really "see" the ends of the sample. In a macroscopic 
sample in which the (single-particle) relaxation length, 
determined, e.g., by the inelastic scattering length, does 
not exceed the size of the system, finite-sizes effects, cor- 
responding to a length scale smaller than the former, 
are naturally smeared out. In the following sections we 
consider nanowire samples that have a nano-meter scale 
cross section, with its circumference sufficiently smaller 
than the relaxation length. Here, we concentrate on the 
cylindrical geometry, imposing additionally a rotational 
(cylindrical) symmetry. We also assume that the system 
is extended to infinity, or (by taking only two of four end 
surfaces into account) periodic in the remaining direc- 
tion. A symptom of the effects we discuss in this section 
may be observed experimentally in a transport measure- 
ment analogous to the one in Ref. 31 . 



Spin-to-surface locking on the cylindrical 
surface 



The protected surface state of a topological insulator 
is often cited with another adjective "helical" . The word, 
helical, stems from a specific feature, often referred to as 
spin-to-momentum locking, 32 that the helical state ex- 
hibits in momentum space. Here, we highlight another 
characteristic of the helical surface state, the "spin-to- 
surfacc locking" , which manifests in real space, and when 
the surface is curved. The electronic spin in a helical state 
on such a curved surface is shown to be locked in-plane 
to the local tangent of the surface. 8-11 

The spin-to-surfacc locking can be also regarded as a 
consequence of (spin) Berry phase of n. In the case of 
rotationally symmetric (cylindrical) wire, the orbital an- 
gular momentum along the axis of the wire is quantized 
to be half-odd integers. This half-odd integral quantiza- 
tion gaps out the spectrum of electronic motion along the 
wire. The spin-to-surfacc locking leads, indeed, irrespec- 
tive of the presence of rotational symmetry, to opening 
of the Dirac spectrum. 

To be explicit let us consider the continuum limit of 
Eqs. (1) and (2), or an effective k p Hamiltonian at the 
r-point (fc = 0), 



in the cylindrical coordinates: 



= \/ x 2 + y 2 , (f> = arctan 



y 



(6) 



Note that our TI sample occupies the interior of a cylin- 
der of radius R. As shown in the Appendix C, any sur- 
face solutions | a)) of Eq. (5) can be expressed as a linear 
combination of the two basis solutions, 

|r+»dv = p(r)\T z +)\r+) dv , 

\r-)) dv = p(r)\r z -)\r-) dv , (7) 

where |r z ±) is an eigenstate of t z with the corresponding 
eigenvalue ±1 and 



|r±) dv - 



1 

T2 



± e i4>/ 2 



(8) 



are two real spin eigenstates pointing either to the cen- 
trifugal (+r) or to the centripetal (— r) direction. In Eqs. 
(7), p(r) is the radial part of the surface wave function lo- 
calized in the vicinity of the surface of the cylinder, given 
explicitly in Eq. (C12). In Eqs. (7), (8) the subscript 
"dv" is added to make explicit that these spinors are 
double- valued. In terms of |r±)) dv , the surface solution 
| a)) reads 



I")) = a+(<f>)\r+))dv + a-((f))\r-)) dv . 



(9) 



Here, the explicit form of the coefficients a±(cf>) is deter- 
mined by solving the eigenvalue problem for the following 
surface effective Hamiltonian, 



H, 



surf 



= A 



1 

R 



d 



<?X + PzO-y 



i.e., 



where 



H suii a((j)) = Ea{4>), 



cc{<j>) = 



ot+{4>) 

a-(4>) 



(10) 



(11) 



(12) 



Notice here that thanks to the rotational symmetry with 
respect to the axis of the cylinder the orbital angular 
momentum L z is a good quantum number, which can 
be simultaneously diagonalized with H sur f and p z . In 
the following, we focus on such surface eigenstates of L z , 
which can be represented in terms of a.{<j>) introduced in 
Eqs. (11), (12) as 



a(4>) = a LztPz (<f>) 
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= e lL ^ 
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e(p)l + T x m(p) + ATyO-^p^, 
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_ a-M _ 




_a_(o)_ 



-Hbuik 

where m(p) = mo+ m 2P 2 - Here, we focus on the isotropic 
case: m 2A1 = mi and A^ = A for fi = x,y,z. We also 
assume e(p) = 0, for simplicity. We then consider the 
eigenvalue problem for Eq. (4), i.e., 



(13) 



a± (0) is specified by the orientation of the surface crystal 
momentum specified by p z and p^ = L z /R. The corre- 
sponding eigenenergy E of H SUI f is then specified by p$ 
and p z as 



(5) 



E = E{p^ Pz ) =±A 



(14) 
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The state \a)) thus given, and specified by the a(<f>) 
given in Eq. (9), signifies a simultaneous eigenstate 
of -ffbuikj L z and p z , which may be also represented 
\L z ,p z )). Eq. (9) implies that such a state is an 
equal-weight superposition of the centrifugal and the 
centripetal spin components given in Eqs. (8), since 
|a + (0)| = |a_(0)|. This signifies that when an electron is 
on the surface of the cylinder at an angle <j> in the config- 
uration space, its spin state is constrained onto the local 
tangent of the cylinder at this position (spin-to-surface 
locking). While an electron travels around the cylinder 
in the configuration space, the corresponding spin frame 
also completes a 2ir rotation in the spin space. 



B. Half-integral quantization of the orbital angular 
momentum and the resulting finite-size energy gap 

Let us reconsider the statue of the angle (j) in different 
steps of the formulation. In the original bulk effective 
Hamiltonian (4) the angle <\> purely specifies the position 
of an electron in the configuration space. This is also 
the case in its eigenstate |a)). Therefore, \a)) must be 
single-valued with respect to the 27r-rotation of (j), 



S+27T 



|a». 



(15) 



On contrary, <f) in |r±) dv specifies the direction of real 
SU(2) spin. Therefore, |r±)) dv is double-valued with re- 
spect to the 27r-rotation of (f>, 



|r±» 



dv | (/> — >-0+27T 



|r±» 



dv ■ 



(16) 



In Eq. (9) these two boundary conditions are compatible, 
only if 



ot{4> + 2ir) = —at(<fi), 



(17) 



i.e., the coefficients a±((f>) are also anti-periodic. In the 
light of Eq. (13), this requires, 



i. -44 



(18) 



i.e., the orbital angular momentum L z is quantized to be 
half-odd integers. 

Notice also that the double- valuedness of |r±)) dv is not 
essential for the half- integral quantization of L z . One can 
equally employ the single- valued version of Eq. (8), 
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(19) 



which is related to |r±)) dv by a simple phase factor, 

|r±) sv = e^/>±) dv . (20) 



In this single-valued basis the surface effective Hamilto- 
nian acquires an additional phase factor n, the spin Berry 
phase, as 
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(21) 



Then, if one employs the same representation (13) for 
the coefficients a, L z takes formally integral values, 
L z = 0, ±1,±2, The corresponding eigenenergy 

E = E(p < f,,p z ) can be also written formally in the same 
way as in Eq. (14). But in that case, p^ in the same 
formula must be reinterpreted as 



P<t> 



Lz + 1/2 
R ' 



(22) 



We have so far seen that whether one employs the 
double-valued [Eq. (8)] or the single-valued [Eq. (19)] 
basis, one finds, as expected, the same gapped spectrum 
given by Eq. (14) with either i) p^ = L z /R with half- 
odd L z [Eq. (18)], or ii) p^ given as in Eq. (22) with 
L z = 0, ±1, ±2, • • • . The magnitude of the energy gap is 
given by twice of 



En — E 
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(23) 



This energy gap due to spin-to-surface locking, or even- 
tually to the doubling of the original two Dirac cones 
through "side surfaces" of the cylinder, decays only alge- 
braically as a function of (inversely proportionally to) the 
circumference of the cylinder. This enhanced finite-size 
energy gap is in marked contrast with that of the slab 
due to mixing of the two surface wave functions sitting 
mainly on the opposing sides of the slab and separated 
by the bulk energy gap. 



IV. CASE OF THE RECTANGULAR PRISM 
GEOMETRY 

In the previous section, we have considered an ideal- 
ized case of the cylindrical geometry, to demonstrate how 
spin-to-surface locking leads to opening of the finite-size 
energy gap. With the rotational (cylindrical) symmetry 
hypothesized, the cylindrical geometry was best suited 
for analytic considerations of the surface state. Here, we 
attempt to realize an equivalent situation in numerical 
experiments in terms of the tight-binding simulation. For 
that purpose, we consider rather prism-shaped samples 
whose cross section on the plane normal to the axis of the 
(right) prism is a rectangle rather than a circle. From the 
viewpoint of topology, such a rectangular prism shape is 
a natural implementation 11 of the cylinder- like geometry 
on the cubic lattice. 

In addition to that aspect as a substitute of a cylin- 
der, there is also a more positive reason we focus here 
on this rectangular prism geometry. In the last few sec- 
tions, throughout the comparison of slab and cylinder, 
we have seen that preventing the communication of two 
Dirac cones sitting on the opposing sides of the sample 
helps protecting the gaplessness of Dirac cones. We have 
so far discussed such switching on and off of this commu- 
nication channel by changing the system's (global) ge- 
ometry. Here, in this section, a new element comes into 
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cases 


type of the phase 


parity of N z 


size gap; N z dependence gap opening mechanism 


(a) 


WTI 


even 


N^ 1 (iii) doubling of Dirac cones due to confinement 


(b) 


WTI 


odd 


"0" (exponentially small) (i) mixing of the opposing sides through gapped bulk 


(c) 


STI 


irrelevant 


(N z + Nx)^ 1 (ii) spin-to-surface locking 



TABLE II. Three typical behaviors of the finite-size energy gap in the rectangular-prism shaped samples. 






FIG. 2. Surface wave function in the rectangular prism ge- 
ometry [Eq. (24)]; WTI phase (m 2z /rn n = 0.2) with N z 
even. Upper: the square of the wave function \ip(z,x)\ 2 with 
N z — 20, N x = 20 and k y — is plotted in the (z, x)-plane. 
Spin and orbital indices are summed over. A± = Au = 1. 
Lower: \ip(x, y, z)| 2 is plotted in the 3D (x, y, z)-space. The 
front, upper and right surfaces correspond, respectively, to the 
ones normal to (—1,0,0), (0, 1,0) and (0,0, 1). Fixed bound- 
ary condition (FBC) in the z- and ^-directions. PBC in the 
y-direction. 



play, the weak indices. As mentioned in the Intoruduc- 
tion, the weak indices have the potentiality of excluding a 
gapless Dirac cone from a surface oriented in a particular 
direction, i.e., that of the weak vector, v = (z/i , v^, v^). 

Folded surfaces of the rectangular prism geometry are 
more adapted for implementing a weak vector as a means 
for eradicating the "dangerous" gapless channels from 
the targeted side surfaces. Another characteristic of the 
WTI surface state is that it exhibits even number of Dirac 



FIG. 3. Plots of the surface wave function in the rectangular 
prism geometry analogous to FIG. 2. Case of N z odd (N z = 
19). WTI phase. N x = 20, k y = 0, A ± = A n = 1. FBC in 
the z- and ^-directions. PBC in the y-direction. 



cones. These two features combine to make gaplessness of 
the surface state of a prism-shaped WTI a rather subtle 
issue, which depends intricately on the geometry and on 
the nature of weak indices. Depending on the relative 
orientation between the weak vector and the surfaces of 
the rectangular prism and on the size of the prism, non- 
compatibility of the surface wave function with a specific 
boundary condition imposed by the geometry leads to, 
or not to opening of a finite-size energy gap. 

The system we consider here has a shape of rectangular 
prism extended in the y-direction. We assume that the 
prism is infinitely long, or periodic, without end surfaces. 
Each cross section of the system at fixed y is restricted 
to a rectangular area of size N z x N x in the (z, x)-plane: 



1 < z < N z , l<x<N x . 



(24) 
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FIG. 4. Plots of the surface wave function in the STI case 
(ra2 Z /m 2 || = 0.3); plots similar to FIG. 2 and FIG. 3. Here, 
the surface wave function is extended over all the four facets 
of the prism. 



The system has two surfaces (^-surfaces) at x = 1 and 
x = N x normal to x — (1,0,0) and two others (z- 
surfaces) at z = 1 and z = N z normal to z — (0,0,1). 
We assume translational symmetry in the y-direction; ky 
is a good quantum number. As for the anisotropy of bulk 
topological insulators, we consider the case of mass pa- 
rameters with uniaxial-type anisotropy as given in Eq. 
(3). In the WTI phase, this corresponds to the case of 
stacked 2D TI layers piled up in the z-direction. 

In the following, we will mainly focus on the WTI 
phase with a specific weak vector v = (0, 0, 1) normal 
to the £-surfaces. Then, gapless Dirac cones are com- 
pletely eliminated from these surfaces, at least in the 
limit of infinitely large surfaces. In the prism geome- 
try (24), the wave function of the corresponding surface 
state has a finite amplitude only on i:-surfaces, barely 
penetrates into the i-side. The Dirac cones forced to be 
localized in each of the ^-surfaces are subject to a par- 
ticular boundary condition imposed by this combination 
of the prism geometry and the weak vector. Compatibil- 
ity or non-compatibility of the surface wave function with 
this specific boundary condition leads to an even/odd fea- 
ture with respect to N z (width of the ai-surfaces) of the 
finite-size energy gap in the WTI phase. After reviewing 
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FIG. 5. Even/odd feature in the finite-size energy gap 
(case of WTI). The mass parameters are on the (red) line 
mo/m2i = —1 of the phase diagram (FIG. 1), slightly be- 
low (WTI-A case, upper panel) and above (STI case, lower 
panel) the phase boundary at m 2 \\/m2± = 1/4. The gap is 
plotted as a function N z . In the WTI-A case: m2 Z /m 2 \\ =0.2, 
Eo = Eo(Nz) shows an even/odd feature, and for N z even the 
gap scales as ~ (N z + 1)~ . In the STI case: m2±/m 2 || = 0.3, 
a weak even/odd feature for small N z is washed out as N z 
increases, and the gap scales as ~ (N z + N x ) . N x = 20. 
A i = An = 1. 



three typical situations we encounter in the analysis of 
the size gap in the WTI and STI phases, we describe the 
nature of even/odd feature in the spirit of k ■ p approxi- 
mation. 



A. Even/odd feature in the WTI phase 

The three typical situations we investigate are the cases 

of 

• WTI with N z even [case (a)], 

• WTI with N x odd [case (b)], and 

• STI [case (c)]. 

The three cases arc also listed in Table II. In our model, 
Eqs. (1), (2), and in the geometry employed, the three 
situations can be realized by a small change of param- 
eters. As for the concrete choice of parameters, we use 
here the following double standard. 33 We first use the 
"theoretical values" that varies on the lines indicated in 



FIG. 1 for the demonstration of crossover from type (c) 
to type (a), and from type (c) to type (b) behaviors. We 
believe that use of these theoretical values help under- 
standing the nature of the phenomenon in the light of 
the phase diagram. Then, in the actual computation of 
the size gap, we also use "experimental values" of the 
parameters that are deduced from experimental data for 
Bi 2 Se 3 . 21 . 34 

The three situations can be easily contrasted by the 
shape of the surface wave function. In the WTI phase 
(FIG. 2 and FIG. 3) the amplitude of the surface wave 
function concentrates on the two x surfaces. The weak 
vector v is here pointed in the direction z, expels the 
surface state from the sides normal to z. In the STI 
phase (FIG. 4), on contrary, the surface state is extended 
over all the four surfaces. In these figures the square of 
the total amplitude of the surface wave function, 

4 

\^(z,x)\ 2 = J2\^^)\ 2 , (25) 
i=i 

is plotted at each point on a cross section (the system is 
translationally invariant in the y-direction) . 

Let us focus on more detailed structures of the shape of 
the surface wave function in the WTI phase, and compare 
the cases of N z even (FIG. 2) and N z odd (FIG. 3). On 
the two x surfaces, the wave function shows a regular 
pattern, vanishing practically at every other layer, when 
N z is odd, whereas in FIG. 2 it is concave shaped (case 
of N z even) . 

This even-odd feature appears more clearly in the be- 
havior of the finite-size energy gap (see FIG. 5) On the 
(red) line mo/rn2|i = —1 of the phase diagram (FIG. 1), 
slightly below [WTI-A: (0,100)] and above [STI: (1,000)] 
the phase boundary at TO2_l/too = — 1/4 the gap is plot- 
ted as a function N z (the number of stacking layers). In 
the WTI case: m2 Z /m 2 \\ = 0.2, E = E (N Z ) shows 
an even/odd feature, and for N z even the gap scales 
as ~ (N z + l) -1 . In the STI case: m2 Z /m 2 \\ = 0.3, a 
weak even/odd feature for small N z is washed out as N z 
increases, and the gap scales as ~ (N z + A^) -1 . In a 
sense, depending on the parity of the number of stacked 
layers, the system becomes either trivial (gapped, when 
N z even) or gapless (when N z odd). Physically, this 
even/odd feature stems from the fact that WTI can be 
viewed as stacked layers of 2D quantum spin Hall states 
(here, stacked in the z-direction). 

B. Effective surface k ■ p theory 

A single Dirac cone cannot be confined (cf. Klein tun- 
neling). This applies to the STI phase we have consid- 
ered in Sec. Ill, in which any surface state, instead of 
being terminated at the end of a plane, continues to the 
adjacent ones, covering the entire surface. In the WTI 
phase, typically two Dirac cones appear on its surfaces, 
i.e., there are "valleys." In that case, one can confine 
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FIG. 6. Shape of the surface wave function: tight-binding 
model vs. fc-p-approximation. |?/j(z)| 2 , the squared amplitude 
of the surface state wave function at x — 1 (and in the case 
of k z — 0) is plotted for the case of N z even (N z = 30, blue 
points). A continuous red curve is the prediction of k-p theory 
[cf. Eq. (31)]. As in FIG. 2, the mass parameters are chosen 
to be mo/m 2 || = —1, m2±/m 2 || = 0.2; other parameters are 
also set as in FIG. 2. 

them in a finite area of the surface. Let us sketch explic- 
itly how this is possible. 

A typical situation we focus on below is the case in 
which two side faces of the prism is normal to the weak 
vector v, implying that there is no Dirac cone on these 
surfaces. In such a situation, the wave function of the 
WTI surface state has a finite amplitude only on the 
remaining two surfaces parallel to v, barely penetrates 
into the side normal to v. The key observation here is 
that the latter can be regarded as a "boundary condition" 
for the wave function that lives mainly on the primary 
parallel surfaces. 

Let us consider a simple and concrete example. In 
the WTI-A phase, shown in FIG. 1, only two ^-surfaces 
are compatible with the presence of gapless Dirac cones; 
the remaining i-surfaces are normal to v — (0, 0, 1). We 
consider the reciprocal space of a ^-surface, spanned by 
k y and k z ; here, we tentatively disregard the presence of 
^-surfaces, pretending as if the translational symmetry 
in the z-direction is still present. Then, on this k — 
{k y , fc z )-plane, two Dirac points appear in the spectrum 
at fci = (0,0) and at k 2 = (0, tt). The spectrum of the 
rectangular prism is obtained, in a crude approximation, 
by projecting E — E(k y , k z ) in the (k y , fc 2 )-plane onto the 
fc y -axis. When two Dirac cones are superposed in this 
projection, a more careful treatment on the boundary 
condition at the corner to the ^-surfaces is needed (see 
below). 

The £-plane on which we focus is bounded by the z- 
surfaces. Penetration of a surface state into the z-sides 
is incompatible with the weak vector, v — (0, 0, 1). This 
may be described by a boundary condition on the surface 
wave function ij)(y, z) on the £-side, 

^(y,z = 0) = 0, ^(y,z = N z + 1) = 0. (26) 

In the fc • p approximation, the wave function tp(y, z) can 
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be constructed by superposing contributions from one 
valley surrounding a Dirac point at k\ and from another 
located at k 2 . As our system is translationally invariant 
in the y-direction, tp(y, z) is expressed in the form of 



if>(y, z) = e tk y y X (z 



(27) 



where x{ z ) should be chosen to satisfy the boundary con- 
ditions (26). This is allowed only when the y-components 
of fei and of k 2 are identical as k\ = (fc ,fci) and 
&2 = (ko,k 2 ). This is indeed the case in the WTI-A 
phase, where ko = 0, k\ = and k 2 = tt. The superposi- 
tion yields 



X(z) = e 



— p i(ki+Pi)z _ e «(fe2+P2)z 



(28) 



where p\ and p 2 are small displacements from the cor- 
responding Dirac points. Note that this automatically 
satisfies the boundary condition at z = 0. If xi z ) with 
Pi = p 2 = (i.e., the superposition of the wave func- 
tions just at the two Dirac points) is compatible with the 
other boundary condition at z = N z + 1, the resulting 
wave function has the zero energy eigenvalue at k y = k a , 
resulting in the gapless surface states. This occurs typ- 
ically at N z odd, and in the WTI-A phase with fci = 
and k 2 — n. Contrastingly, if finite displacements (i.e., 
PitP2 7^ 0) are necessary to satisfy the boundary condi- 
tion, a finite size gap inevitably appears. Naturally, the 
latter applies to the case of N z even. These two contrast- 
ing behaviors explain the nature of the even/odd feature 
demonstrated in FIG. 5. 

Let us further quantify the case of N z even. To fulfill 
the requirement of Eq. (26) we set p\ = —p 2 = q. The 
boundary condition at z — N z + 1 is satisfied, if 



It is suggestive to apply the above k ■ p effective theory 
to the case of WTI-B and WTI-C phases, (see FIG. 1). 
These two topologically different WTI phases appear on 
the blue line m /m 2 || = —5 in the phase diagram with 
the phase boundary at m 2 ±/m 2 ^ = —1/4. The crossover 
of the finite-size energy gap at the transition between 
these two WTI phases is precisely in parallel with the 
one between STI and WTI-A phases (on the red line: 
wo/?7i2|| = ~~ 1 i n FIG. 1) we have considered so far. In 
the case of WTI-B and WTI-C phases, The constituent 
surface Dirac cones on the k = (k y , k z ) -plane appear at 
ki = (tt, 0) and at k 2 — (0, tt) in the WTI-B phase, and at 
ki = (tt, 0) and at k 2 — (tt, tt) in the WTI-C phase. Here, 
the relative position of the two Dirac cone is essential. In 
the case of WTI-C phase, one can construct the surface 
wave function (28) compatible with the specific boundary 
condition (26) precisely in parallel with the previous case 
of the WTI-A phase, simply by replacing k = with 
k = tt, leading to the same even/odd feature. Notice 
that the surface Dirac cone in the WTI-C phase appears 
in the spectrum of prism geometry E = E pr i sm (k y ) at 

ky = TT. 

In the case of WTI-B phase, the two Dirac cones at 
k\ = (tt, 0) and at k 2 = (0, tt) are projected onto a dif- 
ferent point on the k y axis, making the previous con- 
struction [Eqs. (27), (28)] impossible. This is, of course, 
consistent with the fact that in the WTI-B phase the 
surface states are not confined to the ^-surfaces. This 
observation, in turn, reveals that the relative orientation 
of the two (even number of) Dirac cones in the WTI is 
indeed imposed by the weak indices. On surfaces paral- 
lel to the weak vector v, they must appear in line in the 
direction of v. 



q = ± 



2(N Z + 1) 



(29) 



and n being an odd integer. The lowest energy solution 
with n = 1 determines the energy gap to be, 



Eq 



A 



2(N z + l) ' 



(30) 



i.e., E scales as (N z + 1) _1 for N z even within the range 
of validity of the k ■ p-approximation. Eq. (30) allows 
for comparing the above simple effective theory with the 
calculated spectrum. This is done in FIG. 5 by plotting 
the energy gap obtained by numerical diagonalization of 
the corresponding tight-binding model against the pos- 
tulated scaling of Eq. (30) . 

A similar comparison can be made for the shape of the 
surface wave function. Plugging Eq. (29) with n = 1 
back into Eq. (28) one finds, 



\x(z)f 



4shr 



N z tt 
2(N Z + 1)' 



(31) 



The shape of this envelop function is to be compared with 
the calculated value of the amplitude of the surface state 
cigenspinor at x = 1, which is shown in FIG. 6. 



C. Effects of disorder 

Let us comment here on the robustness of the surface 
states discussed in the previous subsections against dis- 
order. A motivation for this is that since disorder leads 
generally to repulsion of the energy levels, one naturally 
questions whether the finite-size effects discussed so far 
are still meaningful when the size gap is perturbed by 
the effects of level repulsion by disorder. The effects of 
disorder is taken into account by introducing a random 
potential V(r), which obeys a uniform distribution in 
the period \—W/2,W/2\ at each site r of the cubic lat- 
tice, i.e., a scalar random potential, cx 1 in the real and 
orbital spin space, which is also cite-diagonal: 

V = 5^V(r)l®|r)<r|, V(r) e [-W/2, W/2] (32) 

r 

is added to the tight-binding Hamiltonian (1) represented 
in the real space. In Eq. (32) the summation over r 
should be taken over all the lattice sites on the cubic lat- 
tice, r — (x, y, z) with x = 1, 2, • • • , N x , y = 1, 2, • • • , 



and z = 1,2, 



,N Z . In the actual computation we set 
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FIG. 7. Surface wave function in the presence of disorder. 
Comparison of the WTI and STI cases: m2 Z /m 2 || = 0.2 (up- 
per) vs. m2 Z /m 2 || = 0.3 (lower). Here, the simulation is done 
for a system of size, N x x N y x N z = 10 x 10 x 10; i.e., N z is 
even. 




FIG. 8. Finite-size energy gap in the rectangular-prism ge- 
ometry plotted as a function of the "width" N x . Comparison 
between the WTI (blue points) and STI (red points) regimes 
in the case of prism thickness N z odd (N z — 9). The loga- 
rithm of the energy gap Eo is plotted vs. N x for demonstrat- 
ing that Eo = Eo(N x ) decays exponentially, showing actually 
an exponentially damped oscillation in the WTI phase. The 
corresponding solutions of Eq. (B10) are a pair of complex 
numbers (see main text for details). The model parameters 
employed are also given there. 




W = 1, mo = — 1, A± = j4|| = 1 in units of m 2 \\ (which 
is set to be unity). 

In FIG. 7 plots similar to FIG. 2, FIG. 3 and FIG. 
4 performed in the presence of disorder are shown. In 
the upper panel (WTI case) the surface wave function is 
localized mainly on one facet of the prism. This is con- 
trasting to the clean cases (FIG. 2, FIG. 3) and to the 
STI case (lower), in which the surface state is extended 
over all the four facets of the prism. The stripe-shaped 
structure is also still visible, indicating that the surface 
wave functions of a specific shape discussed in the previ- 
ous subsection possess some robustness against disorder. 



D. STI more gapped than WTI !? 

We finally discuss the ^-dependence of the size gap. 
As shown in Table II, there are three different types of 
behaviors in the iV^-dependence of the size gap, each 
corresponding to the three different gap-opening mech- 



FIG. 9. Size dependence of Eo = Eo(N x ) in the case of N z 
odd (N z =9). A plot similar to FIG. 8 but in the case of 
model parameters, yielding as solutions for p in Eq. (B10), 
two real solutions given in the main text. The data points for 
the WTI and STI cases are shown, respectively, in blue and 
in red. 



anisms we have highlighted in this paper. Here, let us 
focus again (c/. FIG. 5) on the (red) line mo/m2\\ = — 1 
in the phase diagram (FIG. 1) slightly above and below 
the phase boundary at TO2_l/ to o = —1/4, and compare 
the STI: (1,000) and WTI-A: (0,001) phases. In the fol- 
lowing demonstrations (FIG. 8, FIG. 9, FIG. 10), how- 
ever, we use slightly different set of parameters inspired 
by the corresponding material parameters of Bi2Se3, 21 
but focus on the same phase boundary between STI and 
WTI-A. Here, the tight-binding parameters are specially 
adjusted 34 to reproduce the band structure in the vicin- 
ity of Z-point obtained by the first-principle calculation. 
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FIG. 10. Size dependence of E — E (N X ) in the case of N z 
even (N z = 10). The mass and velocity parameters are chosen 
to be the same as in the case of FIG. 8. Here, the vertical 
axis for Eq is in the linear scale. For N z even Eq(N x ) shows 
at most a power-law decay, whether the system is in the STI 
or WTI phase (see Table II). The data points for the WTI 
and STI cases are as before shown, respectively, as blue and 
red filled circles. 



The employed parameters are given explicitly as 

mo = -0.1, m 2z = m 2 j_ = 0.1, m 2 \\ = 1, 

A z = A x = 0.1, A {{ = 0.3. (33) 

Here, the parameters are normalized in units of m 2 \\ — 
2.60 eV. This set of parameters corresponds to the case 
of the STI phase. To achieve a weak phase we modify the 
value of m 2 ± in Eq. (33) as m 2 ± — > 0.01. This indeed 
falls on the WTI-A phase in FIG. 1. The spectrum of 
the strong phase is "gapped" , showing a finite-size energy 
gap due to spin-to-surface locking, which decays only al- 
gebraically, E ~ (N z + A^) -1 7^ 0. In the weak phase, 
and in the case of N z odd considered here, the spectrum 
is "gapless" , decaying exponentially as a function of the 
distance ~ N x between the two ideally gapless patches 
(log -Eg cx — N x , E ~ 0). This is indeed a comparison 
of the cases (b) and (c) in Table II. In FIG. 8, the log- 
arithm of the energy gap Eq is plotted vs. N x taking 
into account such an expected exponential decay in the 
WTI-A phase. But here, a systematic deviation from a 
simple exponential decay can be clearly seen, implying 
that this is rather a damped oscillation. 

As mentioned in Appendix B, the magnitude of the 
finite-size energy gap in the slab is directly related to the 
(complex) penetration depth of the surface wave func- 
tion, or pi t2 given in Eq. (B10). One can indeed verify, 



E (N x )oc p{ 



JVx + l 



N x +1 
Pi 



(34) 



Recall that in the WTI-A phase considered here two 
Dirac cones, one at k\ — (0, 0) and the other at k 2 = 
(0, 7r), are well grounded on the i-surfaces. The corre- 
sponding surface wave functions exhibit different pene- 
tration depths at each Dirac point, which are specified 
by Eq. (B10). The solutions of Eq. (B10) at k x = (0,0) 



P = Pi,2(ki) ^0.826 ±0.238 i, (35) 

while they are given by 

/0 = Pi,2(fc 2 )- 0-843 ± 0.166 i, (36) 

at k 2 = (0, 7r), i.e., in the two cases, they become a pair 
of complex numbers. In the slab, the finite-size energy 
gap is k = (k z ,k x )-ieso\ved; Eq = E (k), simply the 
minimal value of which determines the actual magnitude 
of the finite-size energy gap. In the case of rectangular 
prism, contributions from k = ki and from k — k 2 are 
superposed to cope with the boundary condition. Notice 
also that here the surface wave function at k — k\ and at 
k = k 2 are both oscillatory [Eqs. (35) and (36)]. These 
two features combine to give the oscillatory pattern of 
log E in the WTI case in FIG. 8. In the figure, two "the- 
oretical" curves for log-Efj are shown in solid curves for 
comparison, not showing a quantitative agreement with 
the actual data. The two curves correspond to the finite- 
size energy gap given as in Eq. (34) at k = k\ (green) and 
at k — k 2 (cyan) estimated under the hypothesis that 
the system is slab-shaped. The actual iV^-dependence of 
log Eq is somewhat in between. 

FIG. 9 is a plot similar to FIG. 8, making the same 
comparison of the STI and WTI-A phases for the same 
N z odd case except that the model parameters are 
slightly modified from Eq. (33) . We replace one of the ve- 
locity parameters A 2 \\ with A 2x — 0.7, leaving A 2y = 0.3 
(the same value as before). This replacement makes the 
corresponding solutions of Eq. (B10) two real solutions, 
indicating that the surface wave function exhibits a sim- 
ple exponential decay. In the WTI-A phase, we have 
chosen as before m 2 ± — 0.01. The behavior of log£?o in 
the WTI case is qualitatively different from the previous 
case. At the two Dirac points, k = k\ and fc = k 2 , in the 
WTI phase, the solutions of Eq. (B10) are 

pi(fci) ~ 0.821, p 2 (k x ) ~ 0.587 (37) 

at ki = (0, 0), while they are given by 



pi(k 2 ) ~ 0.905, p 2 (k 2 ) ~ 0.532 



(38) 



at k 2 — (0, 7r). The actual magnitude of the size gap 
is determined by the largest value of pi, 2 , which is 
the value of p\ at A; = k 2 . Indeed, the actual in- 
dependence of log Eq approaches to this scaling behavior 
(Eq oc pi(fc2) JVx , shown in a solid straight line in FIG. 9) 
for large enough N x . 

Through these two examples we can convince ourselves 
that in this configuration imposed by the combination of 
the prism geometry and a specific choice of the weak 
vector, which can be achieved by adjusting the direction 
of crystal growth direction with respect to the prism, the 
strong topological insulator is qualitatively more gapped 
than a weak topological insulator. 

In the last figure, FIG. 10, we make a comparison be- 
tween the cases (a) and (c) in Table II, in contrast to 
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the previous plots, the ones in FIG. 8, FIG. 9. The 
model parameters are the same as in FIG. 8, but here 
the number N z of stacking layers is even (N z = 10). 
In the STI case the size gap shows a power law decay, 
E ~ (N z + Nx)^ 1 7^ 0, due to spin-to-surfacc locking. 
In the WTI-A phase, the size gap implied by Eq. (30) 
does not scale as a function of N x , but given simply by 



A 2 



2{N Z 



—it = 0.1 x — . 
1) 22 



(39) 



In FIG. 10 this value is indicated as a horizontal grid line 
(in blue) . For sufficiently large value of N x the data looks 
almost constant at a value not much far from the one of 
Eq. (39). 

We have seen so far that from the viewpoint of the 
scaling behavior of finite-size energy gap, the statue of 
the strong and weak phases could be reversed. Here, to 
illustrate this feature we have considered only a very rep- 
resentative range of parameters, but the same feature is 
generic to the vicinity of transitions between STI and 
WTI phases with a suitable choice of the surface direc- 
tions and the number of quintuple layers. 



V. CONCLUSIONS 

We have studied the finite-size energy gap in 3D weak 
and strong topological insulators. Employing the stan- 
dard Wilson-Dirac type effective model, we have devel- 
oped both numerical and analytical considerations. It 
has been demonstrated that anisotropy of the model and 
the geometry of the system are among other model pa- 
rameters crucial elements for determining the qualitative 
nature of the finite-size energy gap. The two elements 
manifest in a correlated manner. The weak topological 
insulator (WTI) has a specific property of (i) expelling 
the gapless surface state from surfaces normal to its weak 
vector V (~ weak indices), i.e., no Dirac cone on the sur- 
face normal to v, (ii) but on surfaces parallel to the weak 
vector, it bears two Dirac cones [more Dirac cones than a 
strong topological insulator (STI)]. We have seen in this 
paper through the study of finite-size effects that these 
two, seemingly competing characteristics of the WTI op- 
erate, in fact, in a cooperative way (c./., k • p-description 
of the surface state in the WTI phase; Sec. IV-B). The 
condition of no Dirac cone on the side normal to V im- 
poses the relative orientation of the two Dirac cones on 
the side parallel to v. The weak indices are also much 
related to the anisotropy of the model parameters. To 
encompass different scaling behaviors of the finite-size 
energy gap, we have manipulated the weak indices by 
varying the model parameters, guided by the phase dia- 
gram shown in FIG. 1. 

Spin-to-surface locking is a characteristic feature of the 
topological insulator surface state, operational both in 
the WTI and STI phases, leading also to a finite-size en- 
ergy gap that exhibits a specific power-law decay as a 
function of the system's linear dimension. Clearly, this is 



more relevant than a usual exponential decay associated 
with the overlap of two surface wave functions, e.g., sit- 
ting on the opposing sides of the slab geometry. By its 
nature the finite-size energy gap due to spin-to-surface 
locking is not effective in the slab, but effective in the 
prism-shaped geometry. In the prism-shaped WTI sam- 
ples, the interplay of these three ingredients; the weak 
vector, the spin-to-surface locking and the rectangular- 
prism geometry leads to intricate finite-size effects, de- 
pending on the model parameters. Three different gap 
opening mechanisms pointed out in this paper: (i) mixing 
of the surface wave functions, (ii) spin-to-surfacc locking, 
and (iii) commensurability with the boundary condition, 
are all effective in determining the intricate size depen- 
dence of the energy gap in the rectangular-prism geome- 
try. 
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Appendix A: Topological numbers 

Notice that our model specified by Eqs. (1), (2) has 
inversion symmetry. This allows us to find the strong 
and weak Z 2 -indices with the use of Fu-Kane's formula. 23 
Here, we mention that in the specific case of e(k) = 
(in most of the analyses in this paper we employ this 
condition for mathematical simplicity) one can introduce 
a Z-type winding number A3. The strong index vq is 
related to N 3 as v n — N 3 mod 2. 

In terms of the periodic table 15,35-39 our starting bulk 
effective Hamiltonian (1) falls on the class AIL This class 
of models has the symmetry, O 2 = — 1, C 2 = and 
= 0, where O, C and T§ represent, respectively, the 
time-reversal, particle-hole and chiral symmetries, and 
in this terminology "0" indicates that the system docs 
not possess that type of symmetry. The periodic table 
says that class All models are characterized by Z2-type 
bulk topoloogical invariants in 3D. For the specific case 
of e(k) — in our model, the symmetry of the model is 
upgraded to the class Dili, i.e., O 2 = — 1, C 2 = 1 and 
= 1, where for the specific Hamiltonian, Eq. (1), C 
and F 5 are given by C = o v T y K and T 5 = r y . This 
symmetry class allows for Z-type bulk topological classi- 
fication in 3D, characterized by a Z-type winding number 
A 3 to be defined below. 

To construct the winding number N 3 explicitly, let us 
first represent the bulk Hamiltonian (1), using an explicit 
matrix representations for the orbital Pauli matrices r x 
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FIG. 11. The winding number N 3 [given in Eq. (A4)], eval- 
uated on a horizontal or vertical line in FIG. 1. In the first 
panel, m2±/m2\\ is fixed at the isotropic point (m,2±/m2\\ = 
1), with mo/m 2 || being varied, while in the remaining pan- 
els mo/m 2 || is at fixed mo/m 2 \\ = — 1 (second panel) and at 
^0/^211 = —5 (third panel). The lines are shown in the same 
color in the phase diagram (see FIG. 1). 



and T y as 



bulk 



m(k)-iP fJl {k)a tl 
m(fe) + iP p (fc)op 



. (A1) 

where we have introduced -P^(fc) = A^sink^. Dividing 
the Hamiltonian by (the magnitude of) its own eigenvalue 
E(k), one can also flatten the spectrum of the Hamilto- 
nian as 



H(k) = 



-tfbulk(fc) 

\E(k)\ 



Q(k) 
Q f (fc) 



(A2) 



where E(k) = ±y / m(k) 2 + P Al (fe) 2 , and 



Q(k) 



m(k) - iPfJJe)^ 
\E{k)\ 



(A3) 



Note that the matrix Q defined above is a 2 x 2 §U(2) 
matrix, satisfying Q^Q — 1 and detQ = 1. Then, one 
can introduce an integral winding number A/3, 40-42 char- 
acterizing the mapping of the 3D Brillouin zone onto this 
SU(2) matrix as 



(A4) 



where T M = Q^dk Q. The integration should be done 
over the entire 3D Brillouin zone. We have evaluated 
this winding number numerically over the entire range of 
parameters shown in FIG. 1 to verify that 



Uq = N 3 mod 2 



(A5) 



indeed holds. The explicit values of A^ in the different 
STI and WTI phases are also shown in FIG. 1. The 
same calculated value is also shown continuously in FIG. 
11 as a function of a control parameter, either TOo/to 2 || 
or m2j_/m 2 || on a few specific lines in FIG. 1. 



Appendix B: Penetration of the surface wave 
function in the slab geometry 

To quantify the surface electronic state in the slab ge- 
ometry, let us concentrate on one surface of the slab. 
Also, we choose this flat surface normal to the x- 
direction. To find the wave function which is localized 
in the vicinity of the surface we divide the bulk Hamil- 
tonian (1) into two parts: 



#buik(fc) = H \\{k\\ 
where fen = (k y , k z ) and 



H±(k x ), 



(Bl) 



(fe|| ) = T x m\\ (fc|| ) + Ty(a y A y sin k y + a z A z sin k z ), 

(B2) 

with my (fen) defined as 

my(fe) = m + 2m2x + 2i r n2y(l — cosfc y ) + 2m2 Z (l — cosfc z ), 

(B3) 

and 



H±(k x ) = —2T x m2 X cosk x 



T y (j x A x sin k x 



(B4) 



This and the following procedure is in parallel with the 
case in which we deal with the continuum model, a more 
standard situation in the context of fc • p approximation, 
discussed in Appendix C, but here we solve the lattice 
model directly without taking the continuum limit. 43 ' 44 
Physically the decomposition (Bl) is based on the picture 
that each (y, z)-plane described by if y (fey) is coupled by 
H±(k x ) to the neighboring layers. In the present geom- 
etry, fen = {k y ,k z ) is a good quantum number. Here, 
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we assume that the system is extended in the half space: 
x > 1, and impose a boundary condition: ip(x = 0) =0. 
A surface solution in such a geometry can be constructed 
by composing a linear combination of base solutions of 
the form, ip(x) — p x ipo (\p\ < !)• For such damped (in- 
stead of plane-wave) solutions, Eq. (B4) modifies to 



-2T x m 2 , 



P 



r v a x A, 



P- P~ 
2i 



(B5) 



In the surface energy spectrum E = E(k\\), protected 
gapless Dirac points can appear at either of the four 
TRIM: fc T RiM = (0,0), (tt,0), (0,tt), (it, it). At such 
TRIM of the surface BZ, the hopping terms in i2||(fc||) 
become inert; 



i?ll(fe|| = A;trim) = T x mu (A;trim)- 



(B6) 



This significantly simplifies the derivation of ^(fcii) at 
k\\ = fexRiM- Notice also that Eq. (B2) with (B3) can be 
regarded as a lattice Hamiltonian for a 2D Z 2 TI with an 
effective mass parameter to 2 d = + 2m 2x . ™2D < 
(to2d > 0) corresponds, respectively, to the non-trivial 
[v = 1) vs. trivial [y = 0) phases, where v is the 2D Z 2 
index. A situation described by this couple of equations 
realizes in the limit N x — > 1. 

Let us construct the surface wave function, 

^^11,^=^0^11), (B7) 
explicitly at fc|| = ^trim- At TRIM, Vo(Ai||) satisfies, 

-ffbuikV'o = [r x m\\ (Aitrim) + H±(p)] V>o = 0, (B8) 
i.e., Vo is a zero-energy eigenstate of 

r x H hulk = m\\(k TmM )-m 2x (p+p~ 1 )+T z <7 x ^(p-p- 1 ). 

(B9) 

Similarly to the case of the continuum model (see Ap- 
pendix C), this zero-energy condition is proven to be 
necessary 40 for constructing a surface solution compat- 
ible with the boundary condition at x = in the form of 
Eq. (B12). Clearly, any of the four simultaneous eigen- 
states of t z and a x , ip±± = |t z ±)|ct 2 ±), is an eigenstate 
of the reduced operator (B9). Then the zero-energy con- 
dition can be used, in turn, to determine p as 



m|| ± Jm 2 -4(m 2 2x -A x /4) 
P = 1 a 7^ = Pl,2 



2(m 2x ± A x /2) 



(B10) 



where 



m\\ = 772 1 (Antrim) = m o (^trim) + 2m 2x . (Bll) 

Here, too^trim) represents the magnitude of bulk en- 
ergy gap at k = fexRiM- In Eq. (B10) the meaning of 
two double signs may need some explanation; the one 
in the numerator is arbitrary, each choice correspond- 
ing to pi.2- The one in the denominator represents + 



for xjjQ = and ip , whereas, the same sign repre- 
sents — for -00 = and ip |_. The structure of Eq. 

(B9) with the understanding that t z <j x — ±1 indicates 
that if p satisfies the zero-energy condition, so docs p^ 1 . 
With a suitable choice of tpo, satisfying both \p\\ < 1 and 
\p 2 \ < 1, the surface solution can be constructed as 



1>(*) = (Pi - 02 M>- 



(B12) 



In a separate paper 45 we study in detail various as- 
pects of the finite-size effects in a slab-shaped sample. 
The magnitude of the finite-size energy gap in the slab is 
determined by the overlap of the two surface wave func- 
tions sitting on the opposing sides of the slab. It is, 
therefore, naturally expected that the magnitude of the 
gap (in a slab of width N x ) is essentially determined by 
the penetration depth, or the amplitude of the wave func- 
tion (B12) at the depth of x = N x . Here, in this model 
one can verify that the correlation of theses two quan- 
tities is a bit stronger than this. The magnitude of the 
size energy gap Eq(N x ) is indeed directly proportional to 
\ip(N x )\ as given in Eq. (34). 



Appendix C: Derivation of the effective surface 
Hamiltonian in the cylinder geometry 

To find the surface effective Hamiltonian on the cylin- 
der in the spirit of k ■ p approximation, 7 ' 11 ' 21 one first 
divides the bulk 3D effective Hamiltonian (4) into two 
parts; one perpendicular, the other parallel to the cylin- 
drical surface: 



H = H ± (p r ) + H\\( P<t> ,p z ), 



(CI) 



where H± = H\ P4>=Pz=0l and p r = —id/dr. H± and H\\ 
read explicitly, 



H± = mj_T x + Ap r T y (a x cos <j> + a y sin <p) 

= t x [m± + iAp r r z f ■ tr] , 
H\\ = m\\T x + At v [p^(- s\n4>a x + cos (/)a y 

= m\\T x + AT y p<f,<p ■ er + p z a z 



(C2) 



PzVz 



where 

777 _L = 777 — 7772 

and 777|| = m 2 (p0 + p 2 z ), with 
.1 d 



dr 2 



ld_ 
r dr 



dz 



(C3) 



(C4) 



(C5) 



We h dso introduced r = (cos <p, sin 0), and <fi = 
(— sin0, cos <f>). 
We then consider a solution of the eigenvalue equation, 



(C6) 
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of the form, <~ e K ( r -- R ) ; i. e ., we set p r = —in (k > 0) 
in Eq. (C2). E± is the value of energy eigenvalue at 
the Dirac point. In order to cope with the boundary 
condition \tp±) r =R = on the surface of the cylinder, 
one can verify that this must be zero (E± = 0). 11,40 This 
implies, 



t x H ± \i/> ± )) =0. 



(C7) 



Notice that in the second line of Eq. (C2) r ■ <x can 
be diagonalizcd by pointing the real-spin spinor in the 
direction of f as Eqs. (8). Then, one can satisfy Eq. 
(C7) by four simultaneous eigenstates of r y and r ■ cr, 
i.e., 

|^)>=p(r)|r a ±)|r±) dv) (C8) 

if k is a solution of 

Ej_ = mj_ ± Ak = m Q - m 2 K 2 ± An = 0. (C9) 

|r±)dv has been given in Eqs. (8). The double sign in Eq. 
(C9) signifies + (— ) when the combination of two signs 
in \T y ±)\f±) in Eq. (C8) are the same (opposite). One 
has to consider a linear combination of the eigenstates of 
the form, 



p{r) 



R) _ e K 2 (r-fl) 



(CIO) 



where K\ and k 2 are solutions of Eq. (C9) with E± = 0, 



i.e., 



±A± \/A 2 + 4m m 2 
4m 2 



«1,2, 



(Cll) 



where the double sign in front of A corresponds to the 
one in Eq. (C9). The second one is arbitrary, each choice 
determining the subscript of k 1-2 - Here, the surface state 
should be localized in the inner vicinity of the surface of 



the cylinder. For that one needs a solution of the form of 
Eq. (C10) with Ki t 2 whose real part both being positive. 
This is in one-to-one correspondence with 

• the choice of + sign in front of A in Eq. (Cll), 
assuming that A/ni2 is positive, and 

• the condition m m 2 < 0. 

Thus, the two basis solutions that span the subspace of 
the surface solutions of Eq. (5) that are also compatible 
with the boundary condition are identified as |r±))d v , 
introduced in Eqs. (7). For preciseness, we normalize 
Eq. (C10) as 



p(r) 



K 2 (r-R) 



nR \Kl — K 2 \ 

Any surface solution |a)) of Eq. (5), satisfying 



(C12) 



ff|||a)) = E\a)), (C13) 
can be expressed as a linear combination of these two 
basis solutions as 



|a)) = a + |r+)) dv + a_|r- )) 



dv • 



(C14) 



or as in Eq. (9). 

Finally, following the prescription of the standard de- 
generate perturbation theory, we consider the secular 
equation for Eq. (C13), i.e., 



((r + |ff|||r+» « r + |if|||r-» 
«r-|ff|||r+» «r-|ff|,|r-» 



a + = E a + 

(C15) 

where we have omitted the subscript "dv" , for simplic- 
ity. We define the coefficient matrix ((r ± |if|||r±)) in 
the secular equation Eq. (C15) as the surface effective 
Hamiltonian H sur f. Noticing the relations such as 

(r±\4>-a\r±)=a y , (C16) 
(f±\a z \f±)=a x , (C17) 

the explicit form of H sur f is found as given in Eq. (10). 
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